/* $Header: p_math.c,v 2.2 01/10/02 10:27:52 przybyls Exp $ */

/***********************************************************************/
/* MODULE NAME:    flowmath                                            */
/* ============                                                        */
/*                                                                     */
/* Author:         Rick                                                */
/*                                                                     */
/* MODULE PURPOSE:  general purpose math routines used for flow/press  */
/*                  calculations (foreground and background)           */
/*                                                                     */
/*                                                                     */
/*                                                                     */
/***********************************************************************/


/*****************************/
/*   GLOBAL INCLUDE FILES    */
/*****************************/


#include <typedefs.h>
#include <m68332.h>
#include <error.h>
#include <p_wkfile.h>
#include <wkfile.h>
#include <p_status.h>
#include <list_mgmt.h>
#include <math.h>
#include <clock_ui.h>   /* import datetime protos */
#include <zone_config.h>
#include <ramps.h>
#include <zone_status.h>
#include <status.h>
#include <keyboard.h>

/*****************************/
/*   LOCAL INCLUDE FILES     */
/*****************************/

#include "pneu_int.h"
#include "proto.h"
#include "p_err.h"

/*****************************/
/*   GLOBAL VARIABLES        */
/*****************************/

#define SMALL_COLUMN   100
#define BIG_IU_SCALE   (16384L)
#define SMALL_IU_SCALE (8L*16384L)

/*****************************/
/*   EXTERNAL VARIABLES      */
/*****************************/

/**begin_proto**/

#pragma SECTION PROG=rom3

/**end_proto**/

 /* required by the floating point routines */

#ifndef mccabe
#pragma ASM

        xdef trap

trap    clr.l   D0
        divu.l   D1,D0

        nop

#pragma END_ASM
#endif

/*****************************/
/*   STATIC VARIABLES        */
/*****************************/

#include "jdp_tab.h"
#include "eta_tab.h"


/***********************************************************************/
/* FUNCTION:                                                           */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:                                                            */
/*                                                                     */
/* RETURNED VALUE:                                                     */
/*                                                                     */
/* GLOBALS USED:                                                       */
/*                                                                     */
/***********************************************************************/
/**begin_proto**/

float GetViscosityExponent( GAS_TYPE gas_type )

/**end_proto**/
{
    return Viscosity_Exp[ (U_INT8) gas_type ];
}

/**begin_proto**/

U_INT16 GetAmbientPressure( U_INT16 Dcm2_Per_IU )

/**end_proto**/
{
     U_INT16 Po;

     switch ( Dcm2_Per_IU )
     {
        case 80:  Po = Atmospheric_Pressure.Dcm2_80;
                  break;

        case 160: Po = Atmospheric_Pressure.Dcm2_160;
                  break;

        case 320: Po = Atmospheric_Pressure.Dcm2_320;
                  break;

        case 640: Po = Atmospheric_Pressure.Dcm2_320/2;
                  break;
              /* default case to make lint happy */
        default:  Po = Atmospheric_Pressure.Dcm2_160;
     }

     return Po;
}

/**begin_proto**/

U_INT16 GetOutletPressure( pCOLUMN_STATUS status )

/**end_proto**/
{
     U_INT16 Po;

     switch ( status->Dcm2_Per_IU )
     {
        case 80:  Po = Atmospheric_Pressure.Dcm2_80;
                  break;

        case 160: Po = Atmospheric_Pressure.Dcm2_160;
                  break;

        case 320: Po = Atmospheric_Pressure.Dcm2_320;
                  break;

        case 640: Po = Atmospheric_Pressure.Dcm2_320/2;
                  break;
              /* default case to make lint happy */
        default:  Po = Atmospheric_Pressure.Dcm2_160;
     }

     return Po;
}

/**begin_proto**/

U_INT16 GetNtpOutletPressure( pCOLUMN_STATUS status )

/**end_proto**/
{
     U_INT16 Po;

     switch ( status->Dcm2_Per_IU )
     {
        case 80:  Po = Standard_Pressure/80;
                  break;

        case 160: Po = Standard_Pressure/160;
                  break;

        case 320: Po = Standard_Pressure/320;
                  break;

        case 640: Po = Standard_Pressure/640;
                  break;
              /* default case to make lint happy */
            default:  Po = Standard_Pressure/160;
     }

     return Po;
}


/**begin_proto**/

U_INT16 CorrectOutletPressure( pCOLUMN_STATUS status, U_INT16 Po )

/**end_proto**/
{

     if ( status->Config->Vacuum_Comp )
     {
          Po = 0;
     }
     else if ( status->Config->Outlet_Pressure_Correction )
     {
          switch ( status->Dcm2_Per_IU )
          {
            case 80:  Po += status->Config->Outlet_Correction;
                      break;

            case 160: Po += status->Config->Outlet_Correction/2;
                      break;

            case 320: Po += status->Config->Outlet_Correction/4;
                      break;

            case 640: Po += status->Config->Outlet_Correction/8;
                      break;

              /* default case to make lint happy */
            default:  Po += status->Config->Outlet_Correction/2;

          }
     }

     return Po;
}


/***********************************************************************/
/* FUNCTION:  CalcInletPresFromJDP                                     */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:        calculate an inlet pressure setpoint                */
/*                                                                     */
/* RETURNED VALUE: a 16bit unsigned int representing inlet pres in IUs */
/*                 this is the Inlet guage pressure!                   */
/*                                                                     */
/* GLOBALS USED:  inlet_pressure_table                                 */
/*                                                                     */
/* CALLS:                                                              */
/*                                                                     */
/* LIMITATIONS:  outlet pres between 8.02 - 15.74 PSI, or 0.           */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcJDP ( U_INT16 Pi, U_INT16 Po_80_dcm2_per_iu )

/**end_proto**/
{

     U_INT16  low_Po, high_Po;
     U_INT16  table_index;

     /* Po indexing and the jdp table is based on a fixed IU scaling */
     /* therefore, Po _must_ be in units of 80 dcm2/iu _always_ !!   */


        /* upper 8 bits of Po becomes the index into */
        /* the two dimensional array of Pi values.   */
        /* later the next 8 bits will be used to     */
        /* interpolate between the two Po entries    */

    table_index = Po_80_dcm2_per_iu >> 8;

        /* since the Po entries do not begin at zero,*/
        /* alter the index to reflect the first entry*/

     table_index = table_index - 0x18;

     low_Po  = tableU ( (U_INT16 *)jdp_table[table_index], Pi );

     table_index++;

     high_Po = tableU ( (U_INT16 *)jdp_table[table_index], Pi );

       /* interpolate between the two Po entries */
     return interpolateU( low_Po, high_Po, (char)((Po_80_dcm2_per_iu) & 0xff) );


}

/***********************************************************************/
/* FUNCTION:  CalcInletPresFromJDP                                     */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:        calculate an inlet pressure setpoint                */
/*                                                                     */
/* RETURNED VALUE: a 16bit unsigned int representing inlet pres in IUs */
/*                 this is the Inlet guage pressure!                   */
/*                                                                     */
/* GLOBALS USED:  inlet_pressure_table                                 */
/*                                                                     */
/* CALLS:                                                              */
/*                                                                     */
/* LIMITATIONS:  outlet pres between 8.02 - 15.74 PSI, or 0.           */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcInletPresFromJDP ( U_INT32 j_delta_p, U_INT16 Po_80_dcm2_per_iu )

/**end_proto**/
{

     U_INT16  low_Po, high_Po;
     U_INT16  table_index;

     /* Po indexing and the jdp table is based on a fixed IU scaling */
     /* therefore, Po _must_ be in units of 80 dcm2/iu _always_ !!   */


        /* upper 8 bits of Po becomes the index into */
        /* the two dimensional array of Pi values.   */
        /* later the next 8 bits will be used to     */
        /* interpolate between the two Po entries    */
    table_index = Po_80_dcm2_per_iu >> 8;

        /* since the Po entries do not begin at zero,*/
        /* alter the index to reflect the first entry*/
     table_index = table_index - 0x18;

     low_Po  = tableU ( (U_INT16 *)inlet_pressure_table[table_index],
                        (U_INT16)j_delta_p );

     table_index++;

     high_Po = tableU ( (U_INT16 *)inlet_pressure_table[table_index],
                        (U_INT16)j_delta_p );

       /* interpolate between the two Po entries */
     return       interpolateU( low_Po,
                                high_Po,
                                (char)((Po_80_dcm2_per_iu) & 0xff) );


}

/***********************************************************************/
/* FUNCTION:                                                           */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:                                                            */
/*                                                                     */
/* RETURNED VALUE:                                                     */
/*                                                                     */
/* GLOBALS USED:                                                       */
/*                                                                     */
/***********************************************************************/
/**begin_proto**/

U_INT16 CalcRealMaxPres( U_INT16 range, U_INT16 Dcm2_Per_IU )

/**end_proto**/
{
    U_INT32 j_delta_p;
    U_INT16 Pi;
    U_INT16 P_outlet;

    switch ( Dcm2_Per_IU )
    {
    case 80 :  j_delta_p = CalcJDP( range/4, Standard_Pressure/80   )*4;
               break;
    default:
    case 160 : j_delta_p = CalcJDP( range/2, Standard_Pressure/80   )*2;
               break;
    case 320 : j_delta_p = CalcJDP( range,   Standard_Pressure/80   );
               break;
    }

    P_outlet = VoltsToAtmPres( 0 );

    switch ( Dcm2_Per_IU )
    {
    case 80:   Pi = (CalcInletPresFromJDP ( j_delta_p/4, P_outlet   ) *4 );
               break;

    case 160:  Pi = (CalcInletPresFromJDP ( j_delta_p/2, P_outlet   ) *2 );
               break;

    case 320:  Pi =  CalcInletPresFromJDP ( j_delta_p,   P_outlet   );
               break;

                    /* default case to make lint happy */
    default:   Pi = (CalcInletPresFromJDP ( j_delta_p/2, P_outlet   ) *2 );

    }

    return Pi;
}

/***********************************************************************/
/* FUNCTION:  ViscosityReferenceTemp                                   */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:   to supply the temperature of the reference viscosity     */
/*                                                                     */
/* RETURNED VALUE:  reference temp in 64ths Kelvin                     */
/*                                                                     */
/* GLOBALS USED:                                                       */
/*                                                                     */
/* CALLS:                                                              */
/*                                                                     */
/* LIMITATIONS:                                                        */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT16 ViscosityReferenceTemp ( GAS_TYPE gas )

/**end_proto**/
{

    switch ( gas )
    {
        case N2:   return N2_Ref_Temp;
        case H2:   return H2_Ref_Temp;
        case He:   return He_Ref_Temp;
        case Argon:   return Ar_Ref_Temp;
        case ArMe: return ArMe_Ref_Temp;
        default:   PneuFatalError(0, BAD_GAS );
    }

    return 1; /* avoid divide by 0 */
}

/***********************************************************************/
/* FUNCTION:  ViscosityReference                                       */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:   to supply a reference viscosity for the current gas      */
/*                                                                     */
/* RETURNED VALUE:  reference viscosity in nPoise                      */
/*                                                                     */
/* GLOBALS USED:                                                       */
/*                                                                     */
/* CALLS:                                                              */
/*                                                                     */
/* LIMITATIONS:                                                        */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT32 ReferenceViscosity ( GAS_TYPE gas)

/**end_proto**/
{

    switch ( gas )
    {
        case N2:      return N2_Ref_Viscosity;
        case H2:      return H2_Ref_Viscosity;
        case He:      return He_Ref_Viscosity;
        case Argon:      return Ar_Ref_Viscosity;
        case ArMe:    return ArMe_Ref_Viscosity;
        default:   PneuFatalError(0, BAD_GAS );
    }

    return 1; /* avoid divide by 0 */
}

/***********************************************************************/
/* FUNCTION:                                                           */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:                                                            */
/*                                                                     */
/* RETURNED VALUE:                                                     */
/*                                                                     */
/* GLOBALS USED:                                                       */
/*                                                                     */
/***********************************************************************/
/**begin_proto**/

void RescaleColumnVelocities( U_INT32 old, U_INT32 new, pCOLUMN_STATUS status )

/**end_proto**/
{
    if ( old == 0 ) return;

    status->Uavg_Setpt = QuadWordMultiplyDivide( status->Uavg_Setpt, new, old );
    status->NTP_Uavg_Setpt = QuadWordMultiplyDivide( status->NTP_Uavg_Setpt, new, old );
    status->Uavg_Actual = QuadWordMultiplyDivide( status->Uavg_Actual, new, old );
    status->NTP_Uavg_Actual = QuadWordMultiplyDivide( status->NTP_Uavg_Actual, new, old );

}
/***********************************************************************/
/* FUNCTION:                                                           */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:                                                            */
/*                                                                     */
/* RETURNED VALUE:                                                     */
/*                                                                     */
/* GLOBALS USED:                                                       */
/*                                                                     */
/***********************************************************************/
/**begin_proto**/

void RescaleColumnFlows( U_INT32 old, U_INT32 new, pCOLUMN_STATUS status )

/**end_proto**/
{
    if ( old == 0 ) return;

    status->Active_Flow = QuadWordMultiplyDivide( status->Active_Flow, new, old );
    status->Desired_Column_Flow = QuadWordMultiplyDivide( status->Desired_Column_Flow, new, old );
    status->Current_V_Flow = QuadWordMultiplyDivide( status->Current_V_Flow, new, old );
    status->NTP_Adjusted_Current_Flow = QuadWordMultiplyDivide( status->NTP_Adjusted_Current_Flow, new, old );

    status->Flow_Setpt = QuadWordMultiplyDivide( status->Flow_Setpt, new, old );
    status->NTP_Flow_Setpt = QuadWordMultiplyDivide( status->NTP_Flow_Setpt, new, old );
    status->Flow_Actual = QuadWordMultiplyDivide( status->Flow_Actual, new, old );
    status->NTP_Flow_Actual = QuadWordMultiplyDivide( status->NTP_Flow_Actual, new, old );

    if ( Initialization_Complete )
       status->Config->Calib.Measured_Flow = QuadWordMultiplyDivide( status->Config->Calib.Measured_Flow, new, old );
}

U_INT8 fPowerOfTwo( float num )
{
    return (U_INT8)( log(num)/log(2) );
}


U_INT16 SlowSqrt( U_INT32 num )
{
     U_INT32   square;
     U_INT32   diff;
     U_INT16   root;
     U_INT16   delta;

     root  = SqrtLU( num );  /* gets within +- 5 */
     delta = 8;

     do
     {
         if ( root * root > num ) root = root - delta;
         else                     root = root + delta;

         delta = delta / 2;
     }
     while ( delta != 0 );

     square = root * root;

     if ( square > num )
     {
         diff = square - num;
         square = ( root - 1 ) * ( root - 1 );
         if ( square > num )
         {
             if ( diff > (square - num) ) root = root -1;
         }
         else
         {
             if ( diff > (num - square) ) root = root -1;
         }
     }
     else
     {
         diff = num - square;
         square = ( root + 1 ) * ( root + 1 );
         if ( square > num )
         {
             if ( diff > (square - num) ) root = root +1;
         }
         else
         {
             if ( diff > (num - square) ) root = root +1;
         }
     }


     return root;
}

/***********************************************************************/
/* FUNCTION:  BuildC1                                                  */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:   perform the math to create a constant value based on     */
/*            pseudo static setpoints that will be used in calculating */
/*            Outlet Velocities.                                       */
/*                                                                     */
/*           2      .7             2                                   */
/*       diam      T      (.0001cm)   10Xdyn/cm2 1e9 nPoise            */
/* C1 = --------------- * ---------   ---------- --------- * IU_Scaling*/
/*        64  len  eta      micron      IU          Poise              */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     X      = IU scaling,  number Pascals per IU                     */
/*                                                                     */
/*   IU_Scaling    = multiplicative offset for integer values,         */
/*              with this correction, upper bits are whole numbers     */
/*                                    lower bits are fractions         */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT32 BuildC1( pCOLUMN_STATUS  status,
                 U_INT16         Dcm2_Per_IU,
                 GAS_TYPE        Gas_Type)

/**end_proto**/
{

     float   opening;
     U_INT32 scaling, old_scaling;
     U_INT8  guess_power;
     float   guess;
     float   diam;
     float   length;

     (void)Gas_Type;

     old_scaling = status->IU_Scaling;

     if ( IsColumnCalibrated( status )  )
     {
         length = (float)status->Config->Calib.Calib_Length;
         diam   = (float)status->Config->Calib.Calib_Diam;
     }
     else
     {
         length = (float)status->Config->Length;
         diam   = (float)status->Config->Diam;
     }

     opening = ((diam - 2 * (float)status->Config->Film_Thickness/COLUMN_FILM_SCALING) * 1e-4 );

     guess =       ( 1 / 64.0 )
                 * ( pow( opening, 2 ) / length )   /* d2/L  in cm     */
                 * Dcm2_Per_IU                      /* for the Pi - Po term */
                 * ( 1e6 * ETA_TABLE_SCALING );     /* vicosity correction (scaled uPoise) */

     guess_power =  fPowerOfTwo( guess );

     if ( guess_power < 21 )  /* room to scale upward */
     {
         guess_power = 21 - guess_power;         /* factors of two we can scale the value by */
         scaling     =  1 << guess_power;        /* grow the scaling */
         guess       = guess * scaling;          /* rescale the value */
     }
     else
     {
         scaling     =  1 << 21;                 /* grow the scaling */
         guess       = guess * scaling;          /* rescale the value */
     }


     status->Vel_IU_Scaling = scaling;

     RescaleColumnVelocities( old_scaling, scaling, status );

     return (U_INT32)guess;

}


/***********************************************************************/
/* FUNCTION:  BuildC2                                                  */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:   perform the math to create a constant value based on     */
/*            pseudo static setpoints that will be used in calculating */
/*            Average Velocities.                                      */
/*                                                                     */
/*              2    .7              2                                 */
/*       3  diam    T       (.0001cm) K  10Xdyn/cm2 1e9 nPoise         */
/* C2 = ----------------  * --------- -- ---------- --------- * IU_Scaling  */
/*      128  len   eta        micron  64   IU        Poise             */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     X      = IU scaling,  number Pascals per IU                     */
/*                                                                     */
/*                                                                     */
/* Note:  also contains a 3/4 multiplicative constant from j.          */
/***********************************************************************/

/**begin_proto**/

U_INT32 BuildC2( pCOLUMN_STATUS  status,
                 U_INT16         Dcm2_Per_IU,
                 GAS_TYPE        Gas_Type)

/**end_proto**/
{

     float opening;
     float   diam;
     float   length;


     (void)Gas_Type;

     if ( IsColumnCalibrated( status )  )
     {
         length = (float)status->Config->Calib.Calib_Length;
         diam   = (float)status->Config->Calib.Calib_Diam;
     }
     else
     {
         length = (float)status->Config->Length;
         diam   = (float)status->Config->Diam;
     }

     opening = ((diam - 2 * (float)status->Config->Film_Thickness/COLUMN_FILM_SCALING) * 1e-4 );


     return  (U_INT32) (
                   ( 3 / 128.0 )                    /*                 */
                 * ( pow( opening, 2 ) / length )   /* d2/L  in cm     */
                 * (float)status->Vel_IU_Scaling

                 * Dcm2_Per_IU                      /* for the Pi - Po term */
                 * ( 1e6 * ETA_TABLE_SCALING ) );   /* vicosity correction (scaled uPoise) */

}


/***********************************************************************/
/* FUNCTION:  BuildC3                                                  */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:   perform the math to create a constant value based on     */
/*            pseudo static setpoints that will be used in calculating */
/*            Volumetric Flow Rates.                                   */
/*                                                                     */
/*                 4 .7           4                                    */
/*     Pntp PI diam T    (.0001cm) K  10Xdyn/cm2 1e9 nPoise 60sec      */
/*C3 = ---------------- * -------- -- ---------- --------- -----* IU_Scaling*/
/*     Tntp 256 len eta    micron  64   IU        Poise      min       */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     Pntp   = 1atm std pressure in IU                                */
/*     Tntp   = 25C  std temperature in Kelvin*64                      */
/*     PI     = 3.141529 ....                                          */
/*     X      = IU scaling,  number Pascals per IU                     */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT32 BuildC3( pCOLUMN_STATUS  status,
                 U_INT16         Dcm2_Per_IU,
                 GAS_TYPE        Gas_Type)

/**end_proto**/
{

     double   opening;
     double   guess;
     U_INT32 scaling, old_scaling;
     U_INT8  guess_power;
     double   diam;
     double   length;

     (void)Gas_Type;

     old_scaling = status->IU_Scaling;

     if ( IsColumnCalibrated( status )  )
     {
         length = (double)status->Config->Calib.Calib_Length;
         diam   = (double)status->Config->Calib.Calib_Diam;
     }
     else
     {
         length = (double)status->Config->Length;
         diam   = (double)status->Config->Diam;
     }

     opening = ((diam - 2 * (double)status->Config->Film_Thickness/COLUMN_FILM_SCALING) * 1e-4 );


     guess =       ( 60.0 * 3.1415926 / 256.0 )                                    /* sec/min, pi     */
                 * ((double)Normal_Temperature / (double)Standard_Pressure )         /* ntp adjustments */
                 * ( pow( opening, 4 ) / length )                                  /* d4/L  in cm     */

                 * Dcm2_Per_IU * Dcm2_Per_IU        /* for the Pi2 - Po2 term */
                 * ( 1e6 * ETA_T_TABLE_SCALING );   /* vicosity correction (scaled uPoise) */


     /* now find a nice big iu_scaling value to keep resolution */

     guess_power =  fPowerOfTwo( guess );

         /* why 21?  flows are 32 bits.  sometimes we need to convert */
         /* inlet flows back to column (scaled) flows.  since the max */
         /* inlet flow is 1000mL, (11 bits) we must only use 32 - 11  */
         /* hence, the largest column flow scale is 21 bits           */

     if ( guess_power < 21 )  /* room to scale upward */
     {
         guess_power = 21 - guess_power;         /* factors of two we can scale the value by */
         scaling     =  1 << guess_power;        /* grow the scaling */
         guess       = guess * scaling;          /* rescale the value */
     }
     else
     {
         scaling     =  1 << 21;                 /* grow the scaling */
         guess       = guess * scaling;          /* rescale the value */
     }


     status->IU_Scaling = scaling;

     RescaleColumnFlows( old_scaling, scaling, status );

     return (U_INT32)guess;
}


/***********************************************************************/
/* FUNCTION: CalcPiFromUavg                                            */
/* ========                                                            */
/*                                                                     */
/* PURPOSE: given a average velocity value (cm/sec *IU_Scaling), calculate  */
/*          calculate a value for inlet pressure.                      */
/*                                                                     */
/* RETURNED VALUE: inlet pressure value in IU                          */
/*                                                                     */
/*                                                                     */
/*********  Base Formula for Average Velocity **************************/
/*                                                                     */
/*            2  .7            2            2                          */
/*      3 diam  T            Pi + 2PiPo + Po               scaling     */
/* Ua = ----------- * ---- * --------------- * (Pi - Po) * constant    */
/*                      .7     2            2              same as C2  */
/*      128 len eta   Tc     Pi +  PiPo + Po                           */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     Tc     = column temperature  (Kelvin*64)                        */
/*     Pi     = inlet pressure in IU                                   */
/*     Po     = outlet pressure in IU                                  */
/*                                                                     */
/*                                                                     */
/*******  Manipulated to solve for Pi  *********************************/
/*                                                                     */
/*                       .7                                            */
/*              Uavg * Tc                                              */
/*  j(Pi-Po) =  ----------                                             */
/*                  C2                                                 */
/*                                                                     */
/* Note: the 3/4 multiplicative constant for j has been lumped into C2.*/
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcPiFromUavg (U_INT32 u_avg, U_INT16 Oven_Temp, U_INT16 Po, pCOLUMN_STATUS status )

/**end_proto**/
{

#define MAX(x,y) (((x)>(y))?(x):(y))

     U_INT32 j_delta_p;
     U_INT32 Pi;
     U_INT16 P_outlet;
     BIT8    ovf;

     P_outlet = CorrectOutletPressure( status, Po );


     j_delta_p = BoundsCheckedQWMD (
                       u_avg,
                       tableU( (U_INT16 *)Eta_Table [ (U_INT8) status->Gas_Type ], Oven_Temp ),
                       status->C2,
                       0xffff,
                      &ovf );

     (void) ovf;

     if (status->Config->Vacuum_Comp == ON)
     {
         if ( j_delta_p > Po )  Pi = j_delta_p - Po;
         else                   Pi = 0;
     }
     else
     {   /* vacuum comp off */


         switch ( status->Dcm2_Per_IU )
         {
             case 80:   Pi = (CalcInletPresFromJDP ( MAX(j_delta_p/4,1), P_outlet) *4 );
                        break;

             case 160:  Pi = (CalcInletPresFromJDP ( MAX(j_delta_p/2,1), P_outlet*2 ) *2 );
                        break;

             case 320:  Pi =  CalcInletPresFromJDP ( MAX(j_delta_p,1),   P_outlet*4 );
                        break;

             case 640:  Pi = (CalcInletPresFromJDP ( MAX(j_delta_p*2,1), P_outlet*8 ) /2 );
                        break;

                    /* default case to make lint happy */
             default:   Pi = (CalcInletPresFromJDP ( j_delta_p/2, P_outlet*2) *2 );

         }
     }

     if ( Pi > 0xffff ) Pi = 0xffff;

     return Pi;

#undef MAX
}

/***********************************************************************/
/* FUNCTION: CalcPiFromVFlow                                           */
/* ========                                                            */
/*                                                                     */
/* PURPOSE: given a volumetric flow (mL/min *IU_Scaling) value, calculate   */
/*          a value for inlet pressure. flow value is normalized to NTP*/
/*                                                                     */
/* RETURNED VALUE: inlet pressure value in IU                          */
/*                                                                     */
/*                                                                     */
/*******  Base Formula for volumetrics flows  **************************/
/*                                                                     */
/*                       4  .7                                         */
/*           Pntp PI diam  T                 2    2                    */
/* Vflow  = ----------------  * ------ * ( Pi - Po ) * scaling constant*/
/*                                 1.7                 same units as C3*/
/*           Tntp 256 len eta    Tc                                    */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     Pntp   = 1atm std pressure in IU                                */
/*     Tntp   = 25C  std temperature in Kelvin*64                      */
/*     PI     = 3.141529 ....                                          */
/*     Tc     = column temperature  (Kelvin*64)                        */
/*     Pi     = inlet pressure in IU                                   */
/*     Po     = outlet pressure in IU                                  */
/*                                                                     */
/*******  Manipulated formula to solve for Pi  *************************/
/*                                                                     */
/*                   1.7                                               */
/*    2      Vflow Tc       2                                          */
/*  Pi  =   ----------- + Po                                           */
/*              C3                                                     */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcPiFromVFlow (U_INT32 v_flow, U_INT16 Oven_Temp, U_INT16 Po, pCOLUMN_STATUS status )

/**end_proto**/
{

     U_INT32 Pi;
     U_INT32 Pi2_Po2;
     U_INT32 Po_sqrd;
     U_INT16 P_outlet;
     BIT8    ovf;


     P_outlet = CorrectOutletPressure( status, Po );

     Po_sqrd = P_outlet * P_outlet;


     Pi2_Po2 = BoundsCheckedQWMD(
                       v_flow,
                       tableU( (U_INT32 *)Eta_T_Table[ (U_INT8) status->Gas_Type ], Oven_Temp ),
                       status->C3,
                       0xffffffff - Po_sqrd,
                      &ovf );

     if ( !ovf )
     {
         Pi = SlowSqrt ( Pi2_Po2 + Po_sqrd );

         if ( Pi <= Po )
         {
              Pi = 0;
         }
         else
         {
              Pi -= Po;
         }
     }
     else
     {
         Pi = 0x10000;
     }

     if ( Pi > 0xffff ) Pi = 0xffff;

     return Pi;

}


/***********************************************************************/
/* FUNCTION: CalcUavg                                                  */
/* ========                                                            */
/*                                                                     */
/* PURPOSE: calculate the average velocity                             */
/*                                                                     */
/* RETURNED VALUE: average velocity (cm/sec)*IU_Scaling                */
/*                                                                     */
/*                                                                     */
/*********  Base Formula for Average Velocity **************************/
/*                                                                     */
/*            2  .7            2            2                          */
/*      3 diam  T            Pi + 2PiPo + Po               scaling     */
/* Uo = ----------- * ---- * --------------- * (Pi - Po) * constant    */
/*                      .7     2            2              same as C2  */
/*      128 len eta   Tc     Pi +  PiPo + Po                           */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     Tc     = column temperature  (Kelvin*64)                        */
/*     Pi     = inlet pressure in IU                                   */
/*     Po     = outlet pressure in IU                                  */
/*                                                                     */
/*                                                                     */
/*******  Whats really done ********************************************/
/*                                                                     */
/*                  2            2                                     */
/*                Pi + 2PiPo + Po                                      */
/*  j_delta_p =   ---------------  ( Pi - Po )                         */
/*                  2            2                                     */
/*                Pi +  PiPo + Po                                      */
/*                                                                     */
/*                                                                     */
/*          C2  j_delta_p                                              */
/*  Uavg =  -------------                                              */
/*                  .7                                                 */
/*                Tc                                                   */
/*                                                                     */
/* Note: the 3/4 multiplicative constant for j has been lumped into C2.*/
/***********************************************************************/
/*                                                                     */
/* LIMITATIONS:   Pi >= 287.5 with  Po = 16.65                         */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT32 CalcUavg ( pCOLUMN_STATUS status,
                   U_INT16        Oven_Temp,
                   U_INT16        P_inlet,
                   U_INT32        Po  )

/**end_proto**/
{

#define MAX(x,y) (((x)>(y))?(x):(y))

     U_INT32 j_delta_p;
     U_INT32 P_outlet;
     U_INT32 Pi;

     /* Zero inlet pressure means velocity = 0 */
     if ( P_inlet == 0 ) return 0;


     P_outlet = CorrectOutletPressure( status, Po );

     Pi = (U_INT32)P_inlet;

     if ( status->Config->Vacuum_Comp )
     {
         j_delta_p = Pi + Po;
     }
     else
     {
          switch ( status->Dcm2_Per_IU )
          {
              case 80 :  j_delta_p = CalcJDP( MAX(Pi/4,1), P_outlet   )*4;
                         break;
              default:
              case 160 : j_delta_p = CalcJDP( MAX(Pi/2,1), P_outlet*2 )*2;
                         break;
              case 320 : j_delta_p = CalcJDP( MAX(Pi,1),   P_outlet*4 );
                         break;
          }
     }

     return QuadWordMultiplyDivide (
                       j_delta_p,
                       status->C2,
                       tableU( (U_INT16 *)Eta_Table [ (U_INT8) status->Gas_Type ], Oven_Temp ) );

#undef MAX
}



/***********************************************************************/
/* FUNCTION: CalcVFlow                                                 */
/* ========                                                            */
/*                                                                     */
/* PURPOSE:  calculate a value for volumetric flow                     */
/*                                                                     */
/* RETURNED VALUE: volumetric flow (mL/min)*IU_Scaling                 */
/*                                                                     */
/*                                                                     */
/*******  Base Formula for volumetrics flows  **************************/
/*                                                                     */
/*                       4  .7                                         */
/*           Pntp PI diam  T                 2    2                    */
/* Vflow  = ----------------  * ------ * ( Pi - Po ) * scaling constant*/
/*                                 1.7                 same units as C3*/
/*           Tntp 256 len eta    Tc                                    */
/*                                                                     */
/*     diam   = column diameter in microns                             */
/*     len    = column length in centimeters                           */
/*     T      = Temperature (Kelvin*64) for viscosity reference        */
/*     eta    = Reference viscosity in nano-poise                      */
/*     Pntp   = 1atm std pressure in IU                                */
/*     Tntp   = 25C  std temperature in Kelvin*64                      */
/*     PI     = 3.141529 ....                                          */
/*     Tc     = column temperature  (Kelvin*64)                        */
/*     Pi     = inlet pressure in IU                                   */
/*     Po     = outlet pressure in IU                                  */
/*                                                                     */
/*******  Manipulated formula to solve for Pi  *************************/
/*                                                                     */
/*                  2    2                                             */
/*           C3  (Pi - Po )                                            */
/* VFlow =   --------------                                            */
/*                1.7                                                  */
/*              Tc                                                     */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT32 CalcVFlow (  pCOLUMN_STATUS status,
                     U_INT16        Oven_Temp,
                     U_INT16        P_inlet,
                     U_INT16        Po  )

/**end_proto**/
{

     U_INT32  Pi;
     U_INT16  P_outlet;

     if ( P_inlet == 0 ) return 0;

     P_outlet = CorrectOutletPressure( status, Po );


     Pi = (U_INT32)P_inlet + Po;

     return   QuadWordMultiplyDivide (
                    (Pi*Pi-P_outlet*P_outlet),
                    status->C3,
                    tableU( (U_INT32 *)Eta_T_Table[ (U_INT8) status->Gas_Type ], Oven_Temp ) );

}


/***********************************************************************/
/* FUNCTION: CalcPiFromPrecWVC                                         */
/* ========                                                            */
/*                                                                     */
/* PURPOSE: calculate a constant flow value for inlet pressure         */
/*                  using vaccuum compensation.                        */
/*                                                                     */
/* RETURNED VALUE: inlet pressure in IU                                */
/*                                                                     */
/*                                                                     */
/*******  Base Formula for Constant Flow *******************************/
/*                                                                     */
/*   V(T1) = V(T2)  Flow is constant at various Temps.                 */
/*                                                                     */
/*                       4  .7                                         */
/*           Pntp PI diam  T                  2    2                   */
/* V(T1)  = ----------------  * ------ * ( Pi1 - Po )                  */
/*                                 1.7                                 */
/*           Tntp 256 len eta    T1                                    */
/*                                                                     */
/*                                                                     */
/*                   2    2            2    2                          */
/*   V(T1) = x * (Pi1 - Po ) = x * (Pi2 - Po ) = V(T2)                 */
/*           --------------        -----------                         */
/*                  1.7                 1.7                            */
/*                T1                  T2                               */
/*                                                                     */
/*                     - or -                                          */
/*                                                                     */
/*                                                                     */
/*                   1.7      2                                        */
/*       2         T2     Pi1                                          */
/*    Pi2       = ---- * ----                                          */
/*                   1.7                                               */
/*                 T1                                                  */
/*                                                                     */
/*                                                                     */
/*     T1     = Temperature (Kelvin*64) for Pi1                        */
/*     T2     = Temperature (Kelvin*64) for Pi2                        */
/*     Pi1    = inlet pressure in IU at T1                             */
/*     Pi2    = inlet pressure in IU at T2                             */
/*     Po     = 0, a vaccuum                                           */
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcPiFromPrevWVC ( pCOLUMN_STATUS status, U_INT16 Oven_Temp, U_INT16 P_25C )

/**end_proto**/
{

       U_INT32 Pi_Setpt_sqrd;
       U_INT16 Pi_Setpt;
       U_INT32 Po;

       Po = GetNtpOutletPressure(status);

       Pi_Setpt_sqrd = QuadWordMultiplyDivide ( ((U_INT32)P_25C + Po ) * (P_25C + Po ),
                                               tableU( (U_INT16 *)Const_Flow_P_Ratio_Table[ (U_INT8) status->Gas_Type ],
                                                                  Oven_Temp ),
                                               CF_P_RATIO_SCALING );

       Po = GetOutletPressure(status);

       Pi_Setpt = SlowSqrt ( Pi_Setpt_sqrd ) ;

       if ( Pi_Setpt <= Po )
       {
           return 0;
       }
       else
       {
           return ( Pi_Setpt - Po );
       }

}


/***********************************************************************/
/* FUNCTION: CalcPiFromPrecWOVC                                        */
/* ========                                                            */
/*                                                                     */
/* PURPOSE: calculate a constant flow value for inlet pressure         */
/*          without using vaccuum compensation.                        */
/*                                                                     */
/* RETURNED VALUE: inlet pressure in IU                                */
/*                                                                     */
/*                                                                     */
/*******  Base Formula for Constant Flow *******************************/
/*                                                                     */
/*   V(T1) = V(T2)  Flow is constant at various Temps.                 */
/*                                                                     */
/*                                                                     */
/*                       4  .7                                         */
/*           Pntp PI diam  T                  2    2                   */
/* V(T1)  = ----------------  * ------ * ( Pi1 - Po )                  */
/*                                 1.7                                 */
/*           Tntp 256 len eta    T1                                    */
/*                                                                     */
/*                                                                     */
/*                   2    2            2    2                          */
/*   V(T1) = x * (Pi1 - Po ) = x * (Pi2 - Po )  = V(T2)                */
/*           --------------        -----------                         */
/*                  1.7                 1.7                            */
/*                T1                  T2                               */
/*                                                                     */
/*                     - or -                                          */
/*                                                                     */
/*                                                                     */
/*                   1.7      2    2                                   */
/*       2    2    T2    (Pi1 - Po )                                   */
/*    Pi2 - Po  = ---- * ------------                                  */
/*                   1.7                                               */
/*                 T1                                                  */
/*                                                                     */
/*                                                                     */
/*     T1     = Temperature (Kelvin*64) for Pi1                        */
/*     T2     = Temperature (Kelvin*64) for Pi2                        */
/*     Pi1    = inlet pressure in IU at T1                             */
/*     Pi2    = inlet pressure in IU at T2                             */
/*     Po     = outlet pressure in IU                                  */
/*                                                                     */
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcPiFromPrevWOVC ( pCOLUMN_STATUS status, U_INT16 Oven_Temp, U_INT16 P_25C )

/**end_proto**/
{

       U_INT32 Pi_sqrd_minus_Po_sqrd;
       U_INT16 Pi_Setpt;
       U_INT32 Po_sqrd;
       U_INT32 Prev_Pi_sqrd;
       U_INT32 Po;

       Po = GetNtpOutletPressure(status);
       Po_sqrd = Po * Po;

       Prev_Pi_sqrd = (U_INT32)( P_25C + Po) * ( P_25C + Po);

       Pi_sqrd_minus_Po_sqrd = QuadWordMultiplyDivide ( Prev_Pi_sqrd - Po_sqrd,
                                tableU( (U_INT16 *)Const_Flow_P_Ratio_Table[ (U_INT8) status->Gas_Type ], Oven_Temp ),
                                CF_P_RATIO_SCALING );

       Po = GetOutletPressure(status);
       Po_sqrd = Po * Po;

       Pi_Setpt = SlowSqrt ( Pi_sqrd_minus_Po_sqrd + Po_sqrd );

       if ( Pi_Setpt <= Po )
       {
           return 0;
       }
       else
       {
           return ( Pi_Setpt - Po );
       }

}


/**begin_proto**/

U_INT16 CalcPiRefFromPi ( pCOLUMN_STATUS status, U_INT16 Oven_Temp, U_INT16 P_setpt )

/**end_proto**/
{

       U_INT32 Pi_sqrd_minus_Po_sqrd;
       U_INT16 Pi_Setpt;
       U_INT32 Prev_Pi_sqrd;
       U_INT32 Po_sqrd;
       U_INT32 Po;

       Po = GetOutletPressure(status);

       Po_sqrd = Po * Po;
       Prev_Pi_sqrd = (U_INT32)( P_setpt + Po) * ( P_setpt + Po);

       Pi_sqrd_minus_Po_sqrd = QuadWordMultiplyDivide ( Prev_Pi_sqrd - Po_sqrd,
                                CF_P_RATIO_SCALING,
                                tableU( (U_INT16 *)Const_Flow_P_Ratio_Table[ (U_INT8) status->Gas_Type ], Oven_Temp ) );

       Po = GetNtpOutletPressure(status);
       Po_sqrd = Po * Po;

       Pi_Setpt = SlowSqrt ( Pi_sqrd_minus_Po_sqrd + Po_sqrd );

       if ( Pi_Setpt <= Po )
       {
           return 0;
       }
       else
       {
           return ( Pi_Setpt - Po );
       }

}


/***********************************************************************/
/***********************************************************************/

/**begin_proto**/

U_INT16 CalcTempFromVFlow (U_INT32 v_flow, U_INT16 Pm, U_INT16 Po, pCOLUMN_STATUS status )

/**end_proto**/
{
     U_INT32 temp;
     U_INT32 Pm2_Po2;
     U_INT32 Po_sqrd;
     U_INT16 P_outlet;
     BIT8    ovf      = 0;

     (void)ovf;

     P_outlet = CorrectOutletPressure( status, Po );

     Po_sqrd = P_outlet * P_outlet;

     Pm2_Po2 = ( (U_INT32)Pm + Po ) *  ( (U_INT32)Pm + Po ) - Po_sqrd;

     temp = BoundsCheckedQWMD ( Pm2_Po2,
                                status->C3,
                                v_flow,
                                0xffffffff,
                               &ovf  );

     temp =  tableU( (U_INT16 *)Temp_Root_Table[ (U_INT8) status->Gas_Type ], temp / (512 * ETA_T_TABLE_SCALING) );

     return temp;
}



/**begin_proto**/

UI_ERR  CalibDiamFromFlow( pCOLUMN_STATUS status )

/**end_proto**/
{

    /*
     *   f = c3 * ( pi2 - po2 ) / etaT
     *
     *   c3 = f * etaT / (pi2 - po2)
     *
     *   d' = d * pow( f * etaT / c3 * (pi2 - po2), 1/4 )
     */


    UI_ERR error;
    U_INT16 oven_init_temp;
    float ft;
    float new_diam;

    oven_init_temp = Active_Wkfile.oven_temp.init_value;

    ft = (2 * (float)status->Config->Film_Thickness/COLUMN_FILM_SCALING) * 1e-4;

    new_diam = pow( ( pow( IsColumnCalibrated(status ) ? status->Config->Calib.Calib_Diam : status->Config->Diam - ft, 4 ) *
                       (float)status->Config->Calib.Measured_Flow /
                       (float)CalcVFlow( status,
                                         oven_init_temp,
                                         CalcColumnInitP( status ),
                                         GetNtpOutletPressure( status ) ) ),
                    0.25 ) + ft;


    if ((error = SetColumnCalibDiam ( status, (INT32)new_diam )) != OK ) return error;

    UpdateColumnCalibLength( status, status->Config->Length );

    SetDateTime( &status->Config->Calib.Date );

   return OK;

}


/**begin_proto**/

UI_ERR CalibDiamFromTime ( pCOLUMN_STATUS status )

/**end_proto**/
{

     /*
      *   d = d' pow( L * t / Uavg , 1/2 )
      *
      */

    UI_ERR error;
    U_INT16 pi, po;
    float ua, length, time, ft;
    float new_diam;

    pi   = CalcColumnInitP( status );
    po   = GetOutletPressure( status );

    ft   = (2 * (float)status->Config->Film_Thickness / COLUMN_FILM_SCALING );

               /* in   cm/sec */
    ua   = (float)CalcUavg ( status, Active_Wkfile.oven_temp.init_value, pi, po ) / (float)status->Vel_IU_Scaling;

               /* in sec */
    time = (float)status->Config->Calib.Holdup_Time / 1000.0;

    if ( IsColumnCalibrated( status ) )
    {
               /* in cm */
        length = (float)status->Config->Calib.Calib_Length;
               /* in iu */
        new_diam = ( status->Config->Calib.Calib_Diam -ft ) * sqrt( length /( time * ua ));
    }
    else
    {
               /* in cm */
        length = (float)status->Config->Length;
               /* in iu */
        new_diam = status->Config->Diam * sqrt( length /( time * ua ));
    }

             /* now put back the diam used up by film thickness */
    new_diam = new_diam + ft;

    if ((error = SetColumnCalibDiam ( status, (INT32)new_diam )) != OK ) return error;

    UpdateColumnCalibLength( status, status->Config->Length );

    SetDateTime( &status->Config->Calib.Date );

    return OK;
}


/**begin_proto**/

UI_ERR CalibLengthFromFlow ( pCOLUMN_STATUS status )

/**end_proto**/
{

    /*
     *   f = c3 * ( pi2 - po2 ) / etaT
     *
     *   c3 = f * etaT / (pi2 - po2)
     *
     *   L' = L * c3 * (pi2 - po2) / f * etaT
     */


    UI_ERR error;
    U_INT32 new_length;
    U_INT16 oven_init_temp;

    oven_init_temp = Active_Wkfile.oven_temp.init_value;


    new_length = QuadWordMultiplyDivide(  IsColumnCalibrated(status) ? status->Config->Calib.Calib_Length:
                                                                       status->Config->Length,
                                         CalcVFlow( status,
                                                    oven_init_temp,
                                                    CalcColumnInitP( status ),
                                                    GetNtpOutletPressure( status ) ),
                                         status->Config->Calib.Measured_Flow );

    if ((error = SetColumnCalibLength ( status, new_length )) != OK ) return error;

    UpdateColumnCalibDiam( status, status->Config->Diam );

    SetDateTime( &status->Config->Calib.Date );

    return OK;
}


/**begin_proto**/

UI_ERR CalibLengthFromTime ( pCOLUMN_STATUS status )

/**end_proto**/
{


     /*
      *   L = L' pow( L' * t * Uavg , 1/2 )
      *
      */

    UI_ERR error;
    U_INT16 pi, po;
    float ua, time, length;
    U_INT32 new_length;


    pi   = CalcColumnInitP( status );
    po   = GetNtpOutletPressure( status );

               /* in   cm/sec */
    ua   = (float)CalcUavg ( status, Active_Wkfile.oven_temp.init_value, pi, po ) / (float)status->Vel_IU_Scaling;


               /* in sec */
    time = (float)status->Config->Calib.Holdup_Time / 1000.0;

    if ( IsColumnCalibrated( status ) )
    {
               /* in cm */
        length = (float)status->Config->Calib.Calib_Length;
               /* in iu */
        new_length  = (U_INT32)(status->Config->Calib.Calib_Length * sqrt( time * ua / length ));
    }
    else
    {
               /* in cm */
        length = (float)status->Config->Length;
               /* in iu */
        new_length  = (U_INT32)(status->Config->Length * sqrt( time * ua / length ));
    }


    if ((error = SetColumnCalibLength ( status, new_length )) != OK ) return error;

    UpdateColumnCalibDiam( status, status->Config->Diam );

    SetDateTime( &status->Config->Calib.Date );

    return OK;
}


/**begin_proto**/

UI_ERR CalibLengthAndDiam ( pCOLUMN_STATUS status )

/**end_proto**/
{

#define MAX(x,y) (((x)>(y))?x:y)

    float  jdp, Pi2_Po2;
    float  Tc, dcm2;
    float  eta, Pi, Po;
    float  Pntp, Tntp;
    float  tm, fntp, ft;
    UI_ERR   error;
    float    len, diam;

    Tc      = (float)Active_Wkfile.oven_temp.init_value / 64;
    Tntp    = (float)Normal_Temperature / 64;

    dcm2    = status->Dcm2_Per_IU;

    Po      = (float)GetOutletPressure( status ) * dcm2;
    Pi      = (float)CalcColumnInitP( status ) * dcm2 + Po;
    Pntp    = (float)GetNtpOutletPressure ( status ) * dcm2;

    Pi2_Po2 = ( Pi * Pi )  -  ( Po * Po );

    eta     = ((float)tableU(
                          (U_INT16 *)Eta_Table [ (U_INT8) status->Gas_Type ],
                          Active_Wkfile.oven_temp.init_value                   ) / ETA_TABLE_SCALING ) / 1e6;

    tm      = (float)status->Config->Calib.Holdup_Time / 1000;

    fntp    = (float)status->Config->Calib.Measured_Flow /( 60 * status->IU_Scaling); /* mL/sec, !mL/min: cgs */

    ft = (2 * (float)status->Config->Film_Thickness/COLUMN_FILM_SCALING);  /* in microns */

    if ( status->Config->Vacuum_Comp )
    {
        jdp = Pi + Po;  /* abs pres */
    }
    else
    {
         jdp = ( (Pi*Pi) + (2*Pi*Po) + (Po*Po) ) * ( Pi - Po ) / ( (Pi*Pi) + (Pi*Po) + (Po*Po) );
    }

    len = pow ( ( fntp  * 9.0 * Tc * Pntp * tm * tm * jdp * jdp ) /
                ( 64.0 * Tntp * 3.141592 * eta * Pi2_Po2 ), 1.0/3.0 );

    diam = ( pow ( (128.0 * len * len * eta ) / (3 * tm * jdp ), .5 ) ) * 1e4 + ft;

    if ((error = SetColumnCalibLength ( status, (INT32)(len  + 0.5) )) != OK ) return error;
    if ((error = SetColumnCalibDiam   ( status, (INT32)(diam + 0.5) )) != OK ) return error;

    SetDateTime( &status->Config->Calib.Date );

    return OK;

#undef MAX
}


/**begin_proto**/

#pragma SECTION UNDO

/**end_proto**/
